rm(list = ls())
gc()
set.seed(63296)
packages <-c("readstata13", "raster", "stargazer", "ggplot2", "spdep", 
             "tseries", "rgdal", "splm", "reshape2","spatialreg")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR DIRECTORY HERE")


# load community shape file
shape <- sf::st_read("./Datasets/spatial_models_data/municipalities_km.shp") %>%
  mutate(gmina = tolower(gmina)) %>%
  filter(woj %in% c('KA', 'CZ', 'OP', 'BB'))

# treatment
treatment <- read.csv("./Datasets/spatial_models_data/map.csv") %>%
  dplyr::select(locality_volha, strikes, priests_continuous, subbotnik_compliance) 

# change label
treatment$strikes <- treatment$strikes
treatment$sabotage <- treatment$subbotnik_compliance

# turn subbotnik upside down to mean sabotage
treatment$sabotage <- treatment$sabotage*-1

# merge both data sets
shape = shape %>% left_join(treatment, by = c('gmina' = 'locality_volha'))

# throw out treatment data
rm(treatment)



####################################
#### Spatial models for strikes ####
####################################

# Restrict data to complete cases
shape_strikes <- shape %>% filter(!is.na(shape$strikes))

#create adjacency matrix
nb_strikes <- poly2nb(shape_strikes, row.names = shape_strikes$gmina)

# This creates the adjacency matrix; here we use a binary matrix ('style="B"') 
# We could also row-standardize ('style="W"') d rather than a weighted matrix
lw_strikes <- nb2listw(nb_strikes, zero.policy=TRUE, style="B" )

# formulas
f1 <- scale(strikes) ~ priests_continuous + factor(woj)

# basic OLS model
model1_strikes_ols <- lm(f1, data = shape_strikes)

# spatial error model
model1_strikes_sem <- errorsarlm(f1, data=shape_strikes, lw_strikes, tol.solve=1.0e-30, zero.policy = T)

# spatial lag model
model1_strikes_slm <- lagsarlm(f1, data=shape_strikes, lw_strikes, tol.solve=1.0e-30)



#####################################
#### Spatial models for sabotage ####
#####################################

# Impute missing sabotage with mean
shape_sabotage <- shape %>% 
  mutate(sabotage = ifelse(is.na(sabotage), mean(sabotage, na.rm=TRUE), sabotage))

#create adjacency matrix
nb_sabotage <- poly2nb(shape_sabotage, row.names = shape_sabotage$gmina)

# This creates the adjacency matrix; here we use a binary matrix ('style="B"') 
# We could also row-standardize ('style="W"') d rather than a weighted matrix
lw_sabotage <- nb2listw(nb_sabotage, zero.policy=TRUE, style="B" )

# formulas
f1 <- scale(sabotage) ~ priests_continuous + factor(woj)

# basic OLS model
model1_sabotage_ols <- lm(f1, data = shape_sabotage)

# spatial error model
model1_sabotage_sem <- errorsarlm(f1, data=shape_sabotage, lw_sabotage, tol.solve=1.0e-30, zero.policy = T)

# spatial lag model
model1_sabotage_slm <- lagsarlm(f1, data=shape_sabotage, lw_sabotage, tol.solve=1.0e-30)


# output regression
stargazer(model1_strikes_sem, model1_strikes_slm, model1_sabotage_sem, model1_sabotage_slm,
          dep.var.labels = c("Strikes", "Sabotage"),
          style = "qje",
          covariate.labels = c("Corrupted priests"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01),
          omit = c("Constant", "factor"),
          omit.stat = c("rsq", "f", "ser"), 
          omit.table.layout = "n",
          add.lines = list(c("Fixed effects", "Yes", "Yes", "Yes", "Yes"),
                           c("Controls", "No", "No", "No", "No")),
          out = "PUT YOUR FILEPATH HERE")







